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ABSTRACT 

An implicit method for radiative transfer in SPH is described. The diffusion 
approximation is used, and the hydrodynamic calculations are performed by a 
fully three-dimensional SPH code. Instead of the energy equation of state for 
an ideal gas, various energy states and the dissociation of hydrogen molecules 
are considered in the energy calculation for a more realistic temperature and 
pressure determination. In order to test the implicit code, we have performed 
non-isothermal collapse simulations of a centrally condensed cloud, and have 
compared our results with those of finite difference calculations performed by 
MB93. The results produced by the two completely different numerical methods 
agree well with each other. 

Subject headings: hydrodynamics - method: numerical - radiative transfer 
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1. Introduction 

In order to understand the various stages of star formation it is essential to follow 
the exact thermal evolution of a collapsing cloud because some important dynamical pro- 
cesses, for example fragmentation (Silk 1977; Low & Lynden-Bell 1976; Rees 1976; Smith 
1977; Bell & Lin 1994; Lodato & Clarke 2004), are closely related to the thermal evolution. 
The thermal evolution of a collapsing cloud is described by radiation hydrodynamics, and 
inevitably involves a numerical solution. However, it is not easy to implement and run a 
multi-dimensional radiation hydrodynamic code, so various approximate methods have been 
used instead of solving the full radiation hydrodynamic equations directly. 

Yorke (1979, 1980); Masunaga et al. (1998); Masunaga & Inutsuka (2000) used one- 
dimensional hydrodynamic codes including a more realistic approach to radiative transfer. 
Such an approach can track the exact thermal evolution of a collapsing cloud, and is easily 
compared to the observational results. However, in a one-dimensional code it is impossible 
to observe some important dynamics, for example fragmentation. 

On the other hand, a simplified treatment for the radiative transfer is needed in multi- 
dimensional hydrodynamic codes. The Eddington approximation was used by Winkler & 
Newman (1980a,b); Boss (1984); Myhill & Boss (1993), while Larson (1969) used the diffusion 
approximation in his one-dimensional simulations. 

However, the multi-dimensional calculations with radiative transfer published so far have 
been done mostly with codes using the finite difference method. Although smooth particle 
hydrodynamics (hereafter SPH) codes are now used quite commonly in the calculations for 
self-gravitating clouds, so far few attempts 1 to include radiative transfer in an SPH code 
have been published (Lucy 1977; Brookshaw 1985; Whitehouse & Bate 2004). The main 
reason is that the diffusion approximation for treating radiative transfer there has a double 
differential in spatial coordinates which is very difficult to evaluate accurately with SPH 
because the particles occupy in principle arbitrary positions. Brookshaw (1985, hereafter 
B85) derived an equation which bypasses the double differential by converting it to a single 
differential. Cleary & Monaghan (1999, hereafter CM99) extended the B85 method to treat 



1 Whitehouse & Bate (2004) was published after the original work of Viau (2001), so we need to emphasize 
the differences between their approach and ours. First of all, Whitehouse & Bate (2004) concentrated on 
one-dimensional tests while this paper presents a fully three-dimensional calculation. Our test calculations 
are shown in Section 3. Furthermore, a more realistic energy calculation has been used in our code rather 
than the ideal equation of state. Details about the energy calculation are presented in Section 2. On the 
other hand, Whitehouse & Bate (2004) used a more detailed treatment of the radiative transfer which allows 
the radiation and gas temperature to be different. 
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cases with a discontinuous conductivity. 

Another difficulty in the treatment of radiative transfer is the large difference in time 
scales. The radiative time scale is much shorter than the dynamical time scale in a collapsing 
cloud. Therefore, an implicit numerical scheme is needed to treat the thermal and dynamical 
evolutions simultaneously. Viau (2001) developed an effective implicit scheme for radiative 
transfer in a fully three-dimensional SPH code after many unsuccessful attempts. 

We have performed a non-isothermal cloud collapse with the implicit code, and com- 
pared our results with those of Myhill & Boss (1993, hereafter MB93). The comparison 
should be useful because the two methods are completely different from each other. 

The implicit scheme which combines the diffusion approximation and SPH is described 
in section 2 in detail. We have used a more realistic energy calculation instead of an ideal 
equation of state for the determination of temperature, and the detailed procedure for the 
specific internal energy calculation is also given in the same section. The non-isothermal 
cloud collapse simulations and the comparison of our results with those of MB93 are given 
in section 3. The summary is in section 4. 



2. Numerical methods 

2.1. Smoothed Particle Hydrodynamics : SPH 

SPH (Lucy 1977; Gingold & Monaghan 1977) is a grid-free and fully Lagrangian method, 
and therefore has been widely used in gravitational collapse simulations. The Lagrangian 
hydrodynamics with self-gravitation are given by 

§ = — « 

— = --VP-V 2 ®, (2) 
Dt p K J 

m = " V ' v ' (3) 

where D / Dt is the Lagrangian derivative, u is the specific internal energy and other variables 
have their usual meaning. Although there are many alternatives for the SPH formulation of 
equations (1) - (3) (Monaghan 1992; Hernquist & Katz 1989), we have used a common form 
given by 

P = ^rnjWij, (4) 
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z j \Pi Pj j 

Here 11^ is the artificial viscosity given by 

~ a c Sy ijfiij + 

where = h^f^i, h is the smoothing length, v i-7 - = v—v,-, = rj-r,-, c Sjij = (c M +c sJ )/2, 

Py = (Pi + Pi) l^i an d a, /3 and 77 are free parameters 2 . Here c s is the sound speed. Although 
the artificial viscosity has a side-effect in differentially rotating systems (Monaghan 1989; 
Morris & Monaghan 1997; Cha & Whitworth 2003b), it is essential for the treatment of shock 
waves in SPH. W in equations (4) - (6) is a kernel function, and the M4 kernel (Monaghan 
& Lattanzio 1984) has been used in our simulations. For additional details about SPH, the 
reader is referred to the reviews by Benz (1990) and Monaghan (1992). 
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2.2. Diffusion approximation in SPH 

The diffusion approximation has been adopted in our three-dimensional SPH code for 
the treatment of radiative transfer, because it is very hard to trace the exact behavior of 
individual photons in the multi-dimensional hydrodynamic code (but see Oxley & Woolfson 
2003). The diffusion approximation is strictly valid only for regions of great optical depth 
(i.e. optically thick medium). However, the early stage of the collapse of a cloud is optically 
thin because of the very low density and effective cooling. Larson (1969) discussed this point 
in his one-dimensional simulations, and concluded that the diffusion approximation can be 
applied in the isothermal stage of a collapsing cloud, because the diffusion approximation 
essentially keeps the cloud temperature at the boundary value in the early stage of the 
collapse. 

The energy equation with the diffusion term is given by 

Du 1 

— = V-V--V-F, 8 

Dt p p w 



2 We have used a = 1.0, (3 = 2.0 and r\ = 0.1 in all simulations, i] is not a critical free parameter 
(Monaghan 1997) for the artificial viscosity, because normally Vi ^ rj. 
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where F is the radiative flux, and is given by 

4acT 3 
F = 

3K R p 



VT, (9) 



where kr is the Rosseland mean opacity, a is c is the speed of light and ctsb is 

the Stephan-Boltzmann constant. If an effective conductivity, Q(= — 16 3^ p ) is defined, 
equation (8) may be rewritten as 

5H = _^ v .v-iv.(Qvr). (io) 

Equation (10) contains a double derivative, and this V 2 operation is very sensitive to the 
disorder in the particle distribution, such that it may cause numerical noise in the simulation. 
B85 suggested another formulation for the diffusion term to avoid this V 2 operator. We 
describe here briefly the treatment of B85. 

In one-dimension, the diffusion term of equation (10) becomes 

p ( dx dx ~^ ^ dx 2 ) ^ ^ 

The derivation for the first term of equation (11) starts from the Taylor expansions of T(x') 
and Q(x') around T(x) and Q(x), respectively, 

Q(x') = Q( x ) + ( x '- x )^ + ^ x '- x ) 2 ^ + ... : (12) 

TV) = nx) + (x-^ + l(x'-x) 2 ^ + ---- (13) 

From equations (12) and (13) one can derive, 

[Q{x')-Q{x)\[T{x')-T{x)\ = ^-*? d ^% (14) 



and the integral interpolant form of equation (14) becomes 

dQdT = _r IQ(x') - Q(x)][T(x') - T(x)] dWjx - x') ^ 
dx dx J x — x' dx 

For the derivation of the second term of equation (11), we will use equation (13), 



I" T(x') — T(x) dW(x — x') , _ r dTdW(x-x') dx ,_l r d 2 T dW (x - x') ^ 

J x — x' dx J dx dx 2 J dx 2 dx 

(16) 
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and it becomes 

*t = _ 2Q rTw-mm*-*)^ (17) 

dx 2 / x — x' ox 



From equations (15) and (17), the diffusion term of equation (8) in one-dimension becomes 

I (WW + Q fT\ = y m 3 {Qj + Qjm-TfidWq 

p \ dx dx dx 2 J pipj Xi — Xj dxi ' 

and the three-dimensional form becomes 

- (VQ • VT + QV 2 T) = V ^Q* + Qi)^ - T i) {T . _ r ) . ViW (19) 



, PiPj 



With equation (18) or (19), the calculation of V 2 can be avoided, and the energy equation 
becomes finally 



Au 
At 



= I E mj (4 + 4 + n„) (V.-V,) we ^ w ' + » )( V Ii) 

(20) 



This method derived by B85 is valid for a constant or smoothly varing conductivity. 
CM99 suggested a modification to treat media with a discontinuous conductivity. They 
considered the continuity of conductive heat flux across the border of adjacent cells (or 
particles), and found a better expression, g^gf, , for the effective conductivity term instead 
of (Qi + Qj) in equation (20). CM99 also performed various tests with the new formulation, 
and some of them are repeated in the next section. 



2.3. Test of the new formulation: thermal conduction 

Given the approximation involved in the derivation presented above, it is not obvious 
that a solution of equation (20) is also a solution of equation (8). Therefore we have per- 
formed a test of the thermal conduction problem to verify the treatment for the double 
derivative. The thermal diffusion equation is given by 

P^ = V-(*VT), (21) 

where t is time and k is the thermal conductivity. The simplest case for the thermal diffusion 
equation is that k is a constant (i.e. a homogeneous medium). In this case equation (21) 
becomes 

dT k P 

— = — V 2 T. 22 
dt pa 
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Here a simple relation between u and T, u = oT is assumed, and a is the specific heat of 
the medium. The treatments of B85 and CM99 are the same in this simple case, and so 
equation (22) can be rewritten as 

AT 2k y mj Tj - T 3 dW l3 
At pier ^—^ pj Xi — Xj dxi 

in one-dimension. 

For this test, we have used 40 particles in < x < n, and a constant h value for all 
particles, p, a and k are set to unity, and all particles keep their original positions in this 
test, so the situation is very similar to a one-dimensional finite difference simulation. The 
temperature distribution at t = is given by 

T(x,0) = sin a;, (24) 

and the boundary temperature is set to OK. The analytic solution of equation (22) is, 

T(x,t) = e^smx. (25) 

Figure 1 shows the results until t = 2. The solid lines in the figure are the analytic solution 
given by equation (25) at different times, and the dots are the results of our numerical 
calculation. The numerical solution reproduces the analytic solution very well and confirms 
the validity of the B85 formulation. Furthermore we have performed several tests for cases 
with a discontinuous conductivity as presented in CM99. We will show one of them. In 
this test the medium is divided into two parts, and each part has different values for the 
conductivity, so the resultant diffusivity ^= ^ j is different in the left and right part of the 
medium. The initial conditions for this test are shown in Table 1. Refer to section 5.3 of 
CM99 for details of this test. Figure 2 shows the results. The method of B85 as modified by 
CM99 (open circles) is more accurate than the original one by B85 (crosses). We have used 
the the modification suggested by CM99 in our three-dimensional test (See section 3). 



Table 1: Initial conditions for discontinuous conductivity test 



left right 
T T~ 
p 1 1 

K 10 1 

a 1 1 
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Fig. 1. — A numerical solution for the heat conduction equation at t — 0.0, 1.0, 1.5 and 2.0. 
The dots are the numerical solution using equation (23), and the solid lines are the analytic 
solution given by equation (25). There is very good agreement between these two solutions. 
The boundary temperature is set to OK in this simulation. 
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Fig. 2. — Numerical solutions for the heat conduction equation with different conductivity 
at t — 0.0, 0.2, 1.0 and 1.5. The open circles are with the expression of CM99, and the crosses 
are with the original method of B85, respectively. The analytic solution is drawn by a solid 
line. From t = 0.2, the numerical results show small deviations from the analytic solution. 
Generally speaking, the modification of CM99 shows a better agreement than B85. 
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2.4. Implicit scheme for radiative transfer 



The radiative time scale is much shorter than the dynamical time scale in a collapsing 
cloud, so the numerical integration of equation (20) is not easy. An implicit scheme has been 
developed to treat the radiative transfer and the dynamical evolution simultaneously (Viau 
2001), and we explain each step of this implicit scheme. 

1. Define a function T from equation (8), 



where iq is the former step value, and u* is the new value updated by the iteration. 
The second and third terms of equation (26) are known for each time step. The goal 
of the iteration is to find the value of u* (or equivalently T*) which will bring T = 0. 

2. Set the boundary values for temperature, T L and T R . These boundary temperatures 
are initially limited by the temperature range of the opacity table (see section 2.5). 

3. Find ul and ur from T L and T R , respectively. At this stage various energy states and 
the dissociation of hydrogen molecules are considered (see section 2.6). 

4. Find the median value, T* using the bisection method or the Van Wijingaarden- 
Dekker-Brent method (e.g. Press et al. 1992). The convergence speed of the Van 
Wijingaarden-Dekker-Brent method is higher than that of the bisection method. Typ- 
ically ~ 10 iterations are required to solve equation (26) for one particle. This iteration 
is started using u* = u® as the initial guess. Rapid convergence here is a concern since 
we are in a double loop in the code. 

5. Find u* from T* 

6. After all values of u* have been found, we compute T again using equation (26) with 
the new u*. Note that only the first and fourth terms change in equantion (26). 

7. Repeat steps 4-6 until T < tol for all particles, using the latest u* found for each 
particle. 

The tolerance, tol in step 6 is closely related to the resolution (and speed as well) of the 
calculation, and set to 10~ 5 in our simulations. Typically several iterations (< 10) are 
required to reach the desired precision. 




(26) 
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In the procedure explained above, note that only the temperature and specific internal 
energy are updated at every iteration, while the pressure and density are fixed. The Rosse- 
land mean opacity, Kr is also updated by the corresponding temperature and fixed density. 
This whole procedure is executed for each particle separately, i.e. the new temperature is 
determined for each particle while keeping all the other particles at their original tempera- 
ture. The convergence according to the same criterion has to be satisfied for all particles. If 
not, the above procedure is repeated until it is satisfied. 

For comparison, the classical method to solve for the gravitational force for the N- 
body problem involves N 2 /2 steps. Using a tree approach brings this to iVlogiV. Here the 
radiative flux equation has to be solved simultaneously for all particles. When the specific 
energy for one particle is varied, all the other ones are fixed. In summary, this method is in 
N 2 steps, times the number of iterations required. This is the price to pay for the different 
time scales between the radiative and dynamical events. But the method converges well and 
yields very good results. The number of iterations required allows taking larger time steps. 

As an example, we have used a small Altix (Intel Itanium2 processor) for the calcula- 
tions with 50000 particles reported below. The results were obtained with our serial code 
version, and the calculation time is just about an hour. We stopped the calculation when 
the maximum density (this happens usually at the cloud center) reaches 10 5 times the initial 
density. Generally speaking, the RHD calculation takes ~10 times more CPU time than the 
HD calculation without radiative transfer. 

The detailed procedures to derive kr and the specific internal energy from a given 
temperature will be described in sections 2.5 and 2.6, respectively. 



2.5. Rosseland mean opacity 

The Rosseland mean opacity, k r should be determined in order to implement radiative 
transfer in the diffusion approximation. The definition of the Rosseland mean opacity is 
given by 

1 Jo Kv dT au 



Kr dr 



(27) 



where B v is the blackbody function at a frequency v and B is the frequency integrated 
blackbody function (= J °° B u dv). k u in equation (27) is determined by 

Kv = -nQ ext -nr 2 d , (28) 
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where n is the number density, Q ext and are the extinction factor and size of dust grains. In 
the lower temperature range (T < 316K), we have used the model for k v developed by Yorke 
(1979) with the more recent values for the extinction factor, Q ext provided by Preibisch et al. 
(1993). The Rosseland mean opacity for this lower temperature region is shown in Figure 3. 
In the figure, the opacity jumps at T ~ 125K because of the sublimation of the ice mantle 
on dust grains. In the higher temperature region (708K < T < 12500K), we have used the 

10 m -I 



6 o. l 



0.01 



0.001 ^ 1 1 — 1 — 

1 10 100 1000 

T (»K) 

Fig. 3. — The Rosseland mean opacity for IK < T < 708K. For T < 316K, we have used 
the model of Yorke (1979) with the recent values for the extinction factor of Preibisch et al. 
(1993). For 316K < T < 708K, we have used simple interpolated values. The jump around 
T ~ 125K is due to the sublimation of the ice mantles. 

model of Alexander & Ferguson (1994) for kr, which considered the absorption of atomic 
lines (with more than 8 million lines) and molecular lines (with nearly 60 million lines). Grain 
absorption and scattering due to silicates, iron, carbon and SiC have also been considered 
in this model. Alexander & Ferguson (1994) tabulated the opacity using a parameter, R, 
defined by 

R=tL, (29) 

where T| is T/10 6 . Figure 4 shows k,r in the ranges —7 < log R < 1 and 708K < T < 12500K. 
There are two jumps at T ~ 1200K and 2000K due to the sublimation of silicates and 
amorphous carbon, respectively. Therefore, all dust components evaporate at T ~ 2000K, 
and the absorption of molecular and atomic lines becomes dominant above that temperature. 
Each curve in the figure is labeled with the value of log R. For the intermediate temperature 



1000 



5000 10* 

T (°K) 



Fig. 4. — The Rosseland mean opacity as a function of temperature for 708K < T < 12500K 
and —7 < logi? < 1 (Alexander & Ferguson 1994). Each curve is labeled with the value 
of log-R. There are jumps at T ~ 1200K and 2000K. They are due to the sublimation of 
silicates and amorphous carbon, respectively. All dust components evaporate at T ~ 2000K, 
and the absorption of molecular and atomic lines become the dominant coolant. 
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region (316K < T < 708K), we have used simple interpolated values for kr (see Figure 3). 



2.6. Thermodynamics 

Temperature has been used as an independent variable in the iteration explained in 
section 2.4, and the specific internal energy of each particle is evaluated at every step of the 
iteration. The specific internal energy can be determined uniquely for a given temperature, 
density and chemical composition. For a more exact determination, the specific internal 
energies of hydrogen, helium and metals are calculated separately. We also consider the 
dissociation of hydrogen molecules in the energy calculation, but the ionization of hydrogen 
atoms is omitted, because ionization is negligible in the temperature range corresponding to 
our simulations. We now explain how to derive the energy. In this description, X, Y and Z 
denote the mass fractions of hydrogen, helium and metals, respectively. 

The total specific internal energy is given by 

u = u(H) + u(H 2 ) + u(H 2diss ) + u(He) + u(M). (30) 
Here u(H) is the specific internal energy of hydrogen atoms, given by 

u{H) = 6 -Xy^, (31) 
2 ran 

where k is the Boltzmann constant, mu is the mass of the hydrogen atom, and y is the ratio 
of atomic to molecular hydrogen, given by 

where p(H) is the density of atomic hydrogen. In the equilibrium state, y is determined by 



y 2 2.11 



exp 



52490 
T 



(33) 



l-y pX 
(Aller 1963). 

u{H 2 ) is the specific internal energy of molecular hydrogen, given by 

Here E(H 2 ) is the energy of hydrogen molecules, and is composed of three terms, 

E(H 2 ) = ^kT + E rot + E vib . (35) 
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The terms of equation (35) are the translational, rotational and vibrational energies, re- 
spectively. The rotational energy of hydrogen molecules is composed of the contributions of 
ortho— and para— H 2 , and is given by 

hfU rp2 d\nz p , f T 2 dln z 
Er0t ~ k (1 - fo) Zp + foZo ' (36) 

where f Q is the fraction for ortho— H 2 , and is 3/4 in the equilibrium state. It is very hard 
to know the exact ratio between ortho— and para—B. 2 in the star forming core, so we have 
used the equilibrium value in our simulations. z Q and z p in equation (36) are the partition 
functions for ortho— and para— H 2 , respectively, and are given by 

'-j(j + l)0rot~ 



z p = Yl ( 2 i + l ) ex p 



J'=0,2,4,- 
oo 



z = (2j + l)exp 

j=l,3,5,- 



T 



(37) 
(38) 



where 6 rot is 85. 4K. The vibrational energy for hydrogen molecules is given by 

Emh = ex P (6 mb /T)-V (39) 

where 9 vib is 6100K. 

The third term of equation (30) is the dissociation energy of hydrogen molecules, and is 

IXyD 



u(H 2 diss) 



(40) 



2 mu ' 

where D is the energy of dissociation for one molecule, and is 4.4773eV. 

u(He) and u(M) are the energies of helium and metals, respectively, and are given by 



u{He) = \y^, 
y ' 2 Amu 



u{M) 



3 kT 



(41) 



(42) 



2 A m m H ' 

where A m is the mean atomic number of metals, and is set to 16.78 in our simulation 
(Cameron 1968). 

The specific internal energy of the gas is determined uniquely with equations (30) - 
(42) at a given temperature, and vice versa. However, an extra iteration is needed to derive 
the temperature from a given specific internal energy. Therefore, we have used T as an 
independent variable rather than u in our code. 
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3. Test for a nonisothermal collapse 

3.1. Initial conditions 

MB93 developed an FDM code for radiation hydrodynamics, and performed simulations 
for a centrally condensed cloud. They used Cartesian and spherical codes, and the Eddington 
approximation for the treatment of radiative transfer. We have performed the same test of 
MB93 to compare our results. The comparison should be meaningful because the method 
of MB93 and ours are based on two completely different philosophies but deal with exactly 
the same problem. 

The same initial conditions have been used in order to compare the results directly. The 
initial cloud has a mass of l.O87M and a radius of 1.1 x 10 16 cm. Solidbody rotation has 
been imposed around the z-axis, and the angular velocity is 8.2 x 10 _12 s _1 . 

The cloud is initially spherical but centrally condensed, and its density profile is given 

by 

p = i 2 4 2 A r (43) 

a/x 2 + Ay 2 + Az 2 

where pi = 4.28 x 10~ 16 g/cm 3 . To implement this density profile, we have changed the mass 
of the particles according to their position. The mass of an individual particle is given by 

m = =, 44 

^x 2 + Ay 2 + Az 2 

where m 8 = l.O87M /Aftotai, and Aftotai is the total number of particles used in the simulation. 
We have used 50000 particles. If the isothermal collapse stage lasts up to p ~ 10~ 13 g/cm 3 , 
this is a sufficient number of particles 3 to satisfy the numerical Jeans condition (Truelove 
et al. 1997; Bate & Burkert 1997; Truelove et al. 1998; Whitworth 1998; Boss et al. 2000). 
MB93 chose this initial density profile to see the effect of radiative transfer and heating 
immediately, and to mimic the prolateness in star forming cores (Myers et al. 1991). 



3.2. Equation of state 

The temperature of the cloud is set to 10K initially, and the thermal evolution during the 
collapse has been tracked using the implicit radiative transfer method explained in section 
2.4. For comparison, we have performed three tests. The only difference between each test 



3 We assume that all particles have the same mass for checking the numerical Jeans condition. 
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is the calculation of the specific internal energy. Tests 1 and 2 use the energy equation of 
state of an ideal gas for the derivation of the energy, which is given by 



kT 



u 



7-1 fim h ' 

where \i is the mean molecular weight, and is given by 



(45) 



X(l-y) Y 
2 4 



16.78 



(46) 



We have used 7 = 5/3 and 7/5 in Tests 1 and 2, respectively. The composition of the cloud 
is X = 0.70, Y = 0.27, Z = 0.03 and y = 0, so //(= 2.385) is assumed to be constant. In 
Test 3, the more realistic energy calculation explained in section 2.6 has been used. Note 
that \x is not constant in Test 3 due to the variation of y. 



In all tests, the gas pressure, P is derived from the pressure equation of state, 

pkT 



P = 



pm H 



(47) 



3.3. Results 

Figures 5-7 show the results of Tests 1, 2 and 3, respectively. They are snapshots at 
t ~ 0.501t// 4 . In the early stages of the collapse, the cloud contracts isothermally, and an 
elongated core forms in the center immediately due to the initial central condensation. The 
central core starts to trap the radiation and becomes adiabatic. The transition density from 
the isothermal to the adiabatic regimes is ~ 10~ 13 g/cm 3 . This transition density is nearly 
the same in all simulations (See Figure 8). In the adiabatic stage, the temperature of the 
central core tends to increase quickly, while the outer parts of the collapsing cloud is still 
isothermal. The adiabatic core is easily distinguished in the temperature profile in Figures 5 
- 7. The density is flatter in the core than in the outer parts of the cloud, and the infalling 
velocity drops suddenly at the edge of the core to form a shock wave. There is no big 
difference between the results of Tests 1, 2 and 3. However, the temperature of the central 
core in Test 1 is higher than that of Tests 2 and 3. The increase of the central temperature 
in Test 1 is faster than in the other simulations because of the higher 7 value. 

Figure 8 shows the evolution of the central core in the p — T plane. For comparison, 
we have drawn two straight lines which show the slopes for 7 = 5/3 and 7/5. When the 



The freefall time (~ 3.38 x 10 3 yrs) is evaluated with the assumption of a uniform density. 
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Fig. 5.— Results for TEST 1 (7 = 5/3) at t = 0.501%/. The top-left plate shows the 
particle positions near the center of the cloud. An elongated central core can be seen. The 
density profile (top right) of the core is nearly flat, and there is an accretion shock around the 
core in the velocity plot (bottom left). The outer envelope of the cloud remains isothermal 
(~ 10K) in the temperature plot (bottom right). 
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Fig. 6.— The results for TEST 2 (7 = 7/5) at t = 0.5014%. The overall features are very 
similar to those of Test 1, but the temperature of the central core is lower because of the 
smaller 7 value. 
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Fig. 7. — The results for TEST 3 (variable 7) at t — 0.5016t//. They are more similar to 
those of Test 2. 
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central core enters the adiabatic regime, Tests 1 (dots) and 2 (long-dash) show slightly 
different evolutions. The slope of Test 1 is steeper than that of Test 2, so the temperature 
in Test 1 should be higher for the same density. In Test 3 (solid line), the evolution of 
central core follows 7 = 5/3 initially, but it becomes closer to 7 = 7/5 when the density 
becomes greater than ~ 10~ 12 g/cm 3 . This change in slope is due to changes in the energy of 
hydrogen molecules. The rotational energy of hydrogen molecules is not important at very 
low temperatures, so the contribution of the translational energy is dominant. However, as 
the temperature increases, the rotational energy becomes more important, so the 7 value 
gets closer to 7/5. Could one use the p — T relation for the cloud center in Figure 8 and 




10 -15 1Q -14 10 -13 !0-12 10-11 10-10 

LogPc (g/ cm3 ) 



Fig. 8. — Evolution of the density and temperature of the cloud center. Dotted, long- 
dashed and solid lines are the results of Tests 1, 2 and 3, respectively. The temperature of 
the collapsing cloud is T ~ 10K until p c ~ 10~ 13 g/cm 3 , and then increases afterwards. The 
effective 7 value for Test 3 is nearly 5/3 in 10 _13 g/cm 3 < p c < 10 _12 g/cm 3 , and then changes 
to 7/5, due to the excitation of the rotational energy of hydrogen molecules. All tests are 
stopped when the first core starts to expand. 

apply it throughout the cloud, or even to other calculations? Figure 9 shows the p — T 
relation for all particles in Test 3 at t — 0.5309t//. The temperature of some particles in 
the density range of 10~ 14 g/cm 3 ~ 10 _13 g/cm 3 is higher than the boundary value (= 10K). 
Furthermore, there is a temperature dispersion at a given density through the cloud except 
at the lower densities. It is due to the non-spherically symmetric collapse, so particles closer 
to the central core are hotter. If a barotropic equation of state (Bodenheimer 1978; Cha & 
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Whitworth 2003a) were used in the simulation, all particles would lie on a line in the p — T 
plane without any dispersion, because there would be no thermal interaction between the 
hot central core and its surroundings. 



oo 
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Fig. 9. — The p — T relation for all particles in Test 3 at t — 0.5309t//. The overall trend 
coincides well with that of Figure 8, but all particles are not on a line. The dispersion in 
the particle distribution is due to the thermal interaction between the hot central core and 
its surroundings, so some particles in the density range 10~ 14 g/cm 3 ~ 10 _13 g/cm 3 show a 
higher temperature than the boundary value (= 10K). 



3.4. Comparison with MB93 

We compare our results with those of MB93 in Table 2. CC and SC in the table mean 
Cartesian and Spherical codes, respectively. There are small differences in the results, the 
temperatures of SPH are lower than those of MB93. We may presume some reasons for this 
discrepancy. First of all, the treatment for radiative transfer is different. We have used the 
diffusion approximation in the simulations, while MB93 used the Eddington approximation. 
It is not easy to predict the resultant difference due to the different treatments for the 
radiative transfer, but the diffusion approximation may reduce the temperature increase 
(Larson 1969). 

Secondly, there is a difference in the energy calculation of hydrogen molecules. We have 
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used equations (34) - (39), but MB93 used slightly different forms for u(H 2 ). According to the 
u(H 2 ) calculation of Boss (1984) (See Appendix B of Boss (1984)), the transition temperature 
from 7 = 5/3 to 7/5 is 100K. However, in our calculation the transition temperature is 
variable, and ~ 40K in Test 3. Therefore, the temperature increase should be slower in our 
simulation. 



We have presented a fully three-dimensional radiation hydrodynamic code based on 
the SPH method and the diffusion approximation. The difficulty encountered in previous 
attempts of how to treat the double derivative in SPH was solved by using the treatment 
developed by B85 to convert the double derivative to a single one. A thermal conduction 
test shows that this treatment works as expected. 

The second difficulty arises from the large difference between the radiative and dynam- 
ical time scales. The radiative time scale is much shorter than the dynamical time scale 
in a collapsing cloud, especially in the isothermal stage. We have developed a fully three- 
dimensional implicit scheme for dealing with the large difference between the time scales. 

To test our implicit scheme, we performed a nonisothermal cloud collapse of a centrally 
condensed cloud. The same simulation has been performed by MB93 using two FDM codes, 
and we have compared our results with those of MB93. The two numerical methods based 
on two completely different philosophies agree with each other. 

Table 2: The central density and temperature 



4. Summary 



case 



Pc 



= 2.2 x l(r 12 g/cm 3 p c = 1.7 x l(T 12 g/cm 3 



Test 1 
Test 2 
Test 3 
CC 



49.9K 42.6K 

39.2K 34.3K 

41.0K 36.8K 
67.0K 



SC 



56.0K 



Comparison of our results with those of MB93. Here CC and SC mean Cartesian and 
Spherical codes, respectively. 
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